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I.  INTRODUCTION  AND  EXECUTIVE  SUMMARY 


A.  Objective.  Tasks,  and  Phases  j 

( 

f 

The  objective  of  this  project  was  to  study  signal  processing 
methods  for  broadband  direction  finding  and  transient  waveform  analysis 
based  upon  the  technique  of  Vandermonde  matrix  factorization.  There 
were  two  main  tasks. 

Task  1  Study  the  properties  and  behavior  of  the  Vandermonde 
factorization  method  for  finding  the  directions  of  distant  targets  that 
emit  transient  waves. 

Task  2  Investigate  the  use  of  Vandermonde  methods  to  resolve 
modal  components  of  acoustic  transient  signals,  e.g.,  echoes  from 
mechanically  resonant  targets  insonified  by  a  broadband,  active  sonar. 

Research  efforts  to  carry  out  these  tasks  were  begun  in  June  1985 
and  continued  until  October  1986.  Those  efforts  were  organized  into  the 
following  five  phases. 

Phase  1  Formalize  the  theory  of  Vandermonde  factorization,  with 
particular  emphasis  upon  the  specific  applications  named  above. 

Phase  2  Groom  the  mechanics  of  the  Vandermonde  matrix 
factorization  for  matrices  of  the  type  likely  to  arise  in  those 
applications. 

Phase  3  For  the  broadband  direction  finding  application,  generate 
data  sets  to  test  direction  finding  performance  as  a  function  of  physical 
parameters  and  system  order.  (As  further  noted  below,  this  phase  was 
terminated  after  it  was  determined  that  the  array  dimensions,  number  of 
floating  point  operations  required,  and  number  of  weeks  needed  to  carry 
out  the  work  were  too  large  to  be  cost-effective.) 

Phase  4  Test  the  use  of  Vandermonde  factorization  for  modal 
analysis  of  transient  waveforms  of  the  type  seen  in  acoustic  echoes 
from  resonant  objects  excited  by  broadband  active  sonar  transmissions. 


Study  the  effect  of  using  realistic  transmissions  (i.e.,  not  ideal 
impulses)  whose  durations  are  not  infinitesimal  as  compared  with  the 
decay  times  of  the  target  resonances,  and  whose  spectral 
characteristics  might  therefore  interact  with  those  of  the  target. 

Phase  5  Demonstrate  the  application  of  the  Vandermonde 
factorization  to  directional  analysis  of  isolated  broadband  transients 
using  realistic  parameter  sets,  using  the  lowest  order  version  of  the 
broadband  direction  finder. 

The  technical  results  of  each  of  these  phases  of  research  are 
summarized  below  in  Section  II. 

B.  Personnel 

Primary  personnel  engaged  in  the  research  efforts  were  the 
following. 

T.  L.  Henderson  -  Principal  Investigator 

R.  S.  Bailey  -  Research  Engineer 

S.  G.  Lacker  -  Student 

C.  Publication  of  Results 

The  results  of  the  research  effort  were  presented  to  the  scientific 
and  technical  community  as  follows. 

(1)  One  technical  paper  was  published  in  full  in  the  Proceedings  of 
the  Nineteenth  IEEE  Asilomar  Conference  on  Circuits,  Systems,  and 
Computers,  after  having  been  presented  at  that  conference  on  November 
7.  1985. 1 

(2)  One  paper  was  presented  at  the  Fall  1986  meeting  of  the 
Acoustical  Society  of  America.2 

(3)  One  technical  paper  was  prepared  for  submission  to  The  Journal 
of  the  Acoustical  Society  of  America.3 


D.  General  Summaru 


The  Vandermonde  factorization  was  found  to  be  useful  not  only  as  a 
concept  but  as  a  computational  tool  for  applications  involving  spectral 
analysis  of  transients.  Application  to  the  analysis  of  transient 
responses  of  resonant  objects  to  broadband  active  sonar  pulses  (pulse- 
compressed  FM  chirps,  uncompressed  FM  chirps,  and  ideal  bandpass  filter 
transient  responses)  was  investigated,  with  interesting  and,  in  some 
cases,  unexpected  results,  as  summarized  below  in  the  summary  of 
Phase  5. 

The  "signal  subspace"  approach  that  is  implicit  in  the  Vandermonde 
factorization  method  of  transient  analysis  does  have  some  known 
limitations.  One  of  them  is  the  requirement  that  the  signal 
dimensionality  K  (the  number  of  spectral  components  or  "poles")  must  be 
finite-,  indeed.  K  must  be  of  modest  size,  in  theory,  to  permit  the 
necessary  matrix  processing.  Another  drawback  is  one  that  afflicts  all 
of  the  "modern  spectral  analysis"  methods:  The  inherent  nonlinearity  of 
the  analysis  technique  makes  it  difficult  to  predict  performance  against 
signals  for  which  the  method  has  not  been  explicitly  tested,  even  when 
the  new  signals  are  linear  combinations  of  those  for  which  the 
performance  is  already  well  understood.4 

The  broadband  direction  finding  application  was  seen  to  be  useful,  at 
least  for  the  case  of  a  single  source  direction  of  arbitrary  spectral 
output.  Unfortunately,  when  the  number  of  sources  K  becomes  large  one 
must  increase  the  dimension  of  the  signal  processor  correspondingly,  and 
that  can  place  severe  demands  upon  the  aperture  weighting  functions  and 
the  frequency  responses  required  for  signal  conditioning.  However,  for 
multiple  sources  that  emit  isolated,  non-simultaneous  broadband 
transients  (so  that  K  =  1  at  any  given  time  and  a  processor  of  low 
dimensional  order  is  therefore  adequate),  the  results  were  very 


3 


encouraging.  Those  results  were  incorporated  in  one  of  the  technical 
papers  that  was  submitted  for  publication.3 

The  results  of  each  of  the  phases  of  research  are  summarized  below, 
following  a  tutorial  introduction  to  the  basic  theory. 


II.  TUTORIAL  INTRODUCTION 


A.  Vandermonde  Factorization 

Given  a  set  {z^ : k  =  1 K}  of  nonzero,  distinct  numbers  (real  or 

complex)  we  define  the  MxK  Vandermonde  matrix  V  whose  (m,k)th 
element  is  (Z|<)m_1.  This  definition  extends  the  classical  definition  of  a 

square  Vandermonde  matrix,  which  is  well  known  to  be  nonsingular.  It 
follows  that  a  Vandermonde  matrix  can  never  have  deficient  rank  since 
it  contains  a  nonsingular,  square,  Vandermonde  submatrix  having  the 
same  number  of  rows  or  columns  (whichever  is  smaller).  An  MxN  matrix 
A  is  said  to  admit  a  Type  2  Vandermonde  matrix  factorization  of 
order  K  if  A=VC,  where  V  is  an  MxK  Vandermonde  matrix  and  C  is  a 
KxN  constant  matrix.  An  NM  vector  function  a(t)  is  said  to  admit  a 
Type  1  Vandermonde  matrix  factorization  of  order  K  if  a(t)  =  Vc(t)  over 
some  specified  time  domain,  where  c(t)  is  a  Kxi  vector  function. 

B.  Application  to  Transient  Spectral  Analusis 

The  Vandermonde  factorization  is  useful  in  certain  problems  where 
observed  "signals"  (e.g.,  sonar  hydrophone  outputs)  are  known  to  lie 
within  a  signal  subspace  of  finite  dimensionality  which,  when  identified, 
reveals  useful  information  about  the  source(s)  of  the  signals.  For 
example,  suppose  a  transient  waveform  a(t)  admits  a  spectral 
decomposition  of  order  K, 

K 

a(t)  =  £  dk  exp(skt)  ,  (1 ) 

k  =  1 


where  the  unknown  parameters  d^'s  and  s^’s  are  allowed  to  be  complex. 

Then  if  a  data  matrix  A  is  constructed  from  a  set  of  uniformly  spaced 
data  samples  {a(t):  t  =  1 ,2 .  etc.}  as  [A]m  n  =  a(m  +  n-2),  then  A  satisfies 

the  Type  2  factorization  A  =  VC,  where  V  is  the  MxK  Vandermonde  matrix 


generated  from  the  set  {z|<  =exp(sj<) :k  =  1 .2 . K}  and  C  is  the  K*N  matrix 

whose  (k,n)th  element  is  d^Cz^)11"1.  It  follows  that  when  A  is 
constructed  from  sampled  values  of  any  observed  transient  waveform 
a(t)  of  the  form  of  Eq.  (1)  and  its  Vandermonde  factorization  is  then 
determined,  the  set  {z^}  will  be  revealed  in  the  second  row  of  V.  (It  is 
required  that  M>K  in  order  to  achieve  a  unique  factorization  that  reveals 
the  Zfc's.)  The  waveform's  resonant  “poles'*  S|<,  whose  imaginary  parts 

express  the  resonance  frequencies  and  whose  real  parts  specify 
the  exponential  decay  rates,  can  then  be  computed  as  S|<  =  log(Z|<).  The 

amptitude/phase  coefficients  d|<  of  the  individual  modal  resonances  can 
then  be  computed  as  =  [C]^  n  zj<1“n  for  any  selected  value  of  n. 

As  applied  to  the  analysis  of  transient  waveforms,  the  Vandermonde 
factorization  amounts  to  an  extended  formalization  of  classical 
Prony/Hankel  techniques,  and  leads  to  some  useful  insights  regarding 
implementation  (as  exemplified  by  the  Vandermonde  factorization 
procedure  that  is  briefly  described  below  in  the  results  summary  of 
Phase  1). 

C.  Application  to  Broadband  Direction  Finding 


The  broadband  direction  finding  application  is  based  upon  the  use  of 
a  line  hydrophone  that  is  equipped  to  provide  simultaneously  a  set  of  M 

different  outputs  (bm(t)  :m  =  1 ,2 . M},  each  of  which  is  extracted  with  a 

different  aperture  weighting  function  wm(z)  and  a  different  signal 
conditioning  filter  whose  impulse  response  is  hm(t),  i.e., 


where  p(z,t)  denotes  the  incident  acoustic  pressure  at  time  t  and 
position  z  along  the  line  hydrophone's  axis.  To  guarantee  the  desired 
direction  finding  behavior,  the  wm(z)’s  (which  have  to  be  functions  of 
bounded  support  if  the  hydrophone  aperture  is  finite)  and  the  hm(t)'s 
must  obey 


wm.l(z) 

=  d[wm(z)]/dt 

(3) 

-c  hm(t) 

=  d[hm+1(t)]/dt 

(4) 

for  m=1 . M-l ,  where  c  denotes  the  speed  of  sound.  It  follows  that 


bm+i(t) 


hm.,(t)  * 


'  +  00 

{d[wm(z)]/dz}p(z.t)dz 
J  -00 


(5) 


lm+1 


(t) 


+  00 


■00 


{  d[p(z,t)]/ dz }  Wm(z)  dz  .  (6) 


where  integration  by  parts  is  used  in  the  last  step  (capitalizing  upon  the 
fact  that  wm(z)  has  bounded  support).  If  the  incident  sound  consists  of 
a  single  plane  wave  coming  from  a  very  distant  source  whose  direction 
cosine  is  u  (relative  to  the  positive  z  axis),  then  the  acoustic  pressure 
p(z,t)  can  be  expressed  as  a  function  of  the  single  variable  (t  +  uz/c). 
This  means  that 

d[p(z,t)]/dz  =  (u/c)x  e[p(z,t)]/8t  ,  (7) 


which  when  substituted  into  Eq.  (6)  gives  the  result 


bm+1(t)  =  (u/c)hm+1(t)  *  (d/dt) 


p(z,t)  wm(z)  dz 


(8) 


Differentiation  and  convolution  are  commutative,  so  the  time 
differentiation  can  be  applied  to  hm+](t)  instead  of  the  integral.  Using 
Eq.  (4)  and  Eq.  (2),  the  following  simple  result  is  then  obtained: 

bm+1(t)  =  ubmW  W 

for  all  t,  for  m=l,2 . M-1.  This  result  can  be  put  into  vector  form:  Let 

b(t)  denote  the  Mxi  vector  [b i  (t) ,  b2(t) . b|vj(t) ]T  so  that 

b(t)  =  [l.u.  u2 . uM”1  ]T  b,(t)  .  (10) 

It  follows  by  linear  superposition  that  if  there  are  K  different  sound 

sources  with  direction  cosines  U}  ,U2 . u^  .  and  whose  individual 

contributions  to  b -j  (t)  are  denoted  as  c-j  (t) .  C2(t) . C|<(t),  then 

K 

b(t)  =  £  [1 ,  uk,  (uk)2, ....  (u^-1  ]T  ck(t)  ,  (11) 

k=1 

i.e., 

b(t)  =  Vc(t)  ,  (12) 

where  V  is  the  MxK  Vandermonde  matrix  whose  (m,k)th  element  is 

(uk)m_1  and  c(t)  =  [c2(t) . C|<(t)]T.  Therefore,  the  Type  1  Vandermonde 

factorization  of  the  output  vector  b(t)  reveals  the  direction  cosines  of 
all  K  wave  sources  (the  u^'s  appear  in  the  second  row  of  V).  This 
direction  finding  process  works  for  waves  of  arbitrary  spectra  and 


bandwidths.  It  is  required  that  M>K;  i.e.,  M  must  be  at  least  as  large  as 
K*l,  where  K  is  the  number  of  wave  sources  that  are  present. 


111.  RESEARCH  SUMMARIES  FOR  EACH  PHASE 


A.  Phase  1:  Formalization  of  the  Theoru  of  Vandermonde  Factorization 

The  basic  theory  of  Vandermonde  factorization  was  formalized  in 
terms  of  block-Hankel  and  unit-parallelogramic  matrices.  The  details 
are  presented  in  Ref.  1.  The  mathematical  developments  were  quite 
successful,  and  led  to  a  self-contained  set  of  useful  rules  and  formulas 
for  manipulating  the  matrices  and  understanding  the  underlying  structure 
of  the  factorization.  Unfortunately,  time  did  not  permit  the  development 
of  theory  to  predict  performance  in  the  presence  of  noise.  Because  of 
the  inherent  nonlinearity  involved  in  solving  A=VC  for  a  constrained  V 
and  C.  and  the  dependence  of  the  solution  upon  null-eigenvectors  of 
singular  matrices,  for  which  perturbation  theory  is  not  very  well 
developed,  such  an  analysis  would  have  been  a  major  undertaking  beyond 
the  scope  of  the  project. 

B.  Phase  2:  Mechanics  of  Factorization 

1  •  Procedure  for  Type  2  Factorization 

Reference  1  describes  the  procedure  that  we  ultimately  developed 
and  tested  for  obtaining  the  Type  2  Vandermonde  factorization  of 
order  K  of  a  given  MxN  matrix  A.  It  utilizes  the  singular  value 
decomposition  (SVD),  and  is  an  improved  and  extended  version  of  an 
SVD-based  procedure  proposed  in  Ref.  5.  The  procedure  consists  of  the 
following  steps. 

Step  1 :  Solve  the  homogeneous  equation  A^X  =  0  for  an  Mx(M-K) 

matrix  X  whose  columns  are  linearly  independent.  (AH  denotes  the 
conjugate  transpose  of  A.)  In  practice  we  perform  the  SVD  of  X  and 
take  the  singular  vectors  corresponding  to  the  smallest  M-K  singular 
values,  and  use  them  as  the  columns  of  X. 


Step  2:  Find  the  (unique)  unit-parallelogramic  matrix  P  (a  unit- 
parallelogramic  matrix  is  one  whose  main  diagonal  elements  all  have 
unit  value,  and  for  which  all  elements  above  the  main  diagonal  or  below 
the  bottom  diagonal  are  zero)  whose  column  space  matches  that  of  X. 
In  practice  this  is  carried  out  by  performing  the  following  substeps. 

Step  2a:  Perform  forward,  Gaussian  elimination  by  columns  (not  by 
rows  as  is  usually  done)  with  partial  pivoting  (in  accordance  with 
Algorithm  2.12  of  Ref.  6,  modified  for  column-wise  rather  than  row¬ 
wise  elimination),  to  reduce  all  elements  above  the  main  diagonal  to 
zero. 

Step  2b  Perform  backward,  simple  Gaussian  elimination  by 
columns,  without  column  interchange,  to  produce  zeroes  below  the 
bottom  diagonal  without  disturbing  the  zeroes  above  the  main  diagonal. 
The  pivot  elements  are  guaranteed  to  be  nonzero  (see  Appendix  of 
Ref.  5),  as  are  the  main  diagonal  elements  of  the  resulting  matrix,  which 
furthermore  preserves  the  column  space  of  X  by  virtue  of  its  having 
been  created  entirely  by  column  operations. 

Step  2c  Scale  the  columns  of  the  resulting  matrix  to  produce  unit 
values  along  the  main  diagonal.  The  resulting  unit-parallelogramic 
matrix  is  identified  as  P.  This  completes  Step  2. 

Step  3  “Toeplitz- ize"  the  matrix  P  by  averaging  its  element 
values  along  diagonals;  i.e.,  by  replacing  P  by  P,  where  [P]jj  is  the 

average  of  all  elements  [P]mn  such  that  m-n=i-j.  If  the  original 

matrix  A  admits  a  Type  2  factorization  exactly,  i.e.,  if  A  =  VC  with  no 
allowance  for  error  due  to  noisy  data,  then  according  to  Ref.  1  P  will 
already  be  a  Toeplitz  matrix  (in  the  sense  that  the  value  of  its  (i.j)th 
element  depends  only  upon  the  value  of  i-j)  .  In  practice,  however,  P 
will  only  approximate  a  Toeplitz  matrix,  and  the  Toeplitz-izing  step  is 
necessary. 
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Step  4  Define  the  (Cth  order  polynomial 


P(z)  =  1  ♦  pi  z  ♦  P2Z2  +  ...  ♦  p«zK  ,  ( 

whose  coefficients  are  P|  =  [P]^ -| .  Solve  for  its  roots  {z|<:k  =  1,2 . K}. 

Step  5  Construct  V  from  the  set  :  k  =  1 .2 . K}  as 


Mm.k  =  (zk)m‘ 


(14) 


for  m=l,2 . M  and  k=l,2 . K. 


Step  6  Compute  C  =  (VHV)-1  VH  A  . 


The  Type  2  factorization  of  A  is  thus  complete. 


(15) 


Our  procedure  for  performing  a  Type  1  factorization  is  to  transform 
it  into  a  Type  2  factorization  that  has  the  same  Vandermonde  matrix  V, 
using  samples  of  the  data  vector  a(t)  at  t=t1,t2,...,  etc.  In  particular. 
a(t)  =  Vc(t)  implies  A=VC  if  we  either  define 


tAli.j 


[a(t j)]|  5  [C]jj  =  [c(tj)]i 


(16) 


or  define 


J  J 

A  =  Z  a(tj)aH(tj)  i  C  =  I  c(tj)cH(t j) VH  .  0  7) 
j=l  j=l 

Therefore,  to  solve  a  Type  1  factorization  we  construct  the  data 
matrix  A  from  sampled  values  of  the  data  vector  a(t)  in  accordance  with 
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either  Fq.  (16)  or  Eq.  (17),  and  then  compute  the  Type  2  factorization 
A=VC.  This  gives  V,  and  c(t)  can  be  computed  as 


c(t)  =  (VHV)-’  VHa(t)  .  (18) 

3.  General  Summaru  of  Factorization  Experiences 

These  factorization  methods  were  tested  in  a  variety  of  cases  and 
appeared  to  work  quite  well.  However,  neither  the  budget  nor  the  scope 
of  work  of  the  project  supported  the  development  and  publication  of  a 
complete,  self-contained,  well  documented,  and  thoroughly  debugged 
software  package  for  general  distribution.  Instead,  the  method  was 
tested  for  software  parameters  (e.g.,  array  dimensions)  tailored  to  the 
particular  examples  being  investigated,  and  the  procedure  was  conducted 
on  a  step-by-step  basis  to  monitor  performance  and  watch  for 
difficulties  such  as  arithmetic  underflow/overflow.  Standard  software 
libraries  were  used  for  the  SVD  and  polynomial  root-solving  operations 
(which  had  to  accommodate  complex  variables  in  both  the  matrix 
elements  and  the  polynomial  coefficients). 

C.  Phase  3;  Parameter  Study  of  Broadband  Direction  Finding  for  Steadu 
State  Sources 

It  was  intended  to  develop  a  procedure  for  generating  case  study 
data  for  the  broadband  direction  finding  application,  assuming  a  modest 
number  of  wave  sources  that  emitted  broadband,  random  processes  of 
relatively  long  duration,  within  a  noisy  background  of  a  much  larger 
number  of  similar  but  weaker  sources  in  random  directions.  (We  knew 
of  no  other  way  to  generate  realistic  background  noise  for  the  special 
multi-output  hydrophone  required  by  our  construction.)  The  task  turned 
out  to  be  much  more  difficult  than  expected,  due  to  the  proverbial  “curse 
of  dimensionality"  that  sometimes  plagues  signal  processing  efforts. 
The  data  storage  and  file  handling  requirements,  as  well  as  the  number 
of  floating  point  operations  required,  were  finally  determined  to  be 
beyond  the  scope  of  the  project.  Attempts  to  incorporate 


approximations  and  limit  the  observation  interval  were  not  successful 
enough  to  reduce  the  computational  effort  to  an  acceptable  level,  and  the 
effort  was  abandoned.  (Ironically,  an  approximation  technique  that 
might  have  solved  the  problem  was  conceived  shortly  after  the  end  of 
the  project  term.) 


D.  Phase  4:  Studies  of  Transient  Analysis  for  Active  Sonar  with 
Resonant  Targets 

The  Vandermonde  factorization  procedure  was  tried  out  as  a  means 
of  identifying  target  resonances.  In  particular,  a  sonar  target  impulse 
response  of  the  form 

K/2 

h(t)  =  U(t)  x  £  A|(  expf-o^t)  cos(2Ttf|<t-Mp|<)  (19) 

k=1 

was  assumed,  where  U(t)  denotes  the  unit  step.  The  resonance 
frequencies  f|<  and  decay  rates  <X|<  would  be  determined  by  observing  the 
sonar  system  output 

a(t)  =  $(t-r)*p(t)«h(t)  ,  (20) 

where  p(t)  denotes  either  the  transmitted  waveform  (if  a  broadband 
impulsive  pulse  of  extremely  short  duration  is  used)  or  its 
autocorrelation  function  (if  a  time- stretched  pulse  is  transmitted  and 
the  hydrophone  output  is  matched  filtered  for  purposes  of  pulse 
compression)  and  t  denotes  the  two-way  travel  time.  Implicit  in  this 
equation  are  several  simplifying  assumptions: 

(1)  the  target  has  zero  Doppler, 

(2)  the  medium  is  stable  enough  that  it  can  be  modeled  as  a  time- 
invariant  linear  system,  and 


(3)  spreading  loss  and  absorption  factors  are  assumed  to  be 
absorbed  into  the  definition  of  p(t). 

For  analysis  purposes  it  can  be  assumed  that  the  echo  delay,  as 
expressed  by  S(t-r),  will  have  been  removed  in  practice,  at  least 
approximately.  Any  residual  delay  that  is  not  removed  can  be  absorbed 
into  the  definitions  of  A|<  and  <P|<  anyway,  and  does  not  influence  the 
resonance  parameters  fj<  and  that  are  to  be  determined.  If  the 
transmitted  sonar  pulse  has  sufficient  bandwidth  then  p(t)  will 
approximate  a  delta  function,  so  that  a(t)  will  have  the  form  of  Eq.  (19) 
to  a  very  good  approximation;  at  least  it  should  have  the  same  set  of 
resonance  frequencies  f^  and  decay  rates  cfy  ,  perhaps  with  altered 

amplitudes  A^  and  phases  <P|<.  Indeed,  the  error  that  is  incurred  in 
making  this  approximation  can  be  regarded  as  a  form  of  “noise”  in  the 
data  a(t),  and  the  primary  purpose  of  our  study  was  to  determine  the 
robustness  of  the  Vandermonde  factorization  technique  when  applied 
directly  to  the  sonar  output  a(t). 

At  the  risk  of  being  redundant  we  shall  elaborate  further:  Unless 
the  transmitted  pulse  is  extremely  impulsive  (i.e.,  very  intense  and  of 
extremely  short  duration)  to  begin  with,  one  would  like  to  deconvolve 
the  transmitted  sonar  pulse  from  the  sonar  echo  to  recover  an 
unadulterated  version  of  the  target  impulse  response  h(t).  However,  in 
practice  the  transmitted  pulse  must  have  limited  bandwidth  and  true 
deconvolution  is  therefore  impossible.  If  the  active  sonar  system  is 
sophisticated  (and  expensive)  enough,  then  one  can  use  classical  pulse 
compression  as  a  substitute  for  true  deconvolution;  i.e.,  one  can  apply  a 
matched  filter.  The  matched  filter  convolves  the  sonar  echo  with  the 
time-reversed  transmitted  pulse.  Mathematically,  the  effect  is  as  if  the 
autocorrelation  function  p(t)  had  been  transmitted;  hence,  the  assertion 
of  Eq.  (20).  The  failure  of  the  matched  filtering  process  to  perform  a 
true  deconvolution  has  the  result  that  a(t)  is  a  "noisy",  i.e.,  distorted, 
version  of  h(t).  If  the  transmitted  pulse  has  sufficient  bandwidth  to 
encompass  the  target  resonances  adequately,  then  this  distortion  should 


be  small,  and  can  be  modeled  by  a  perturbation  of  the  amplitude  and 
phase  parameters  Ak  and  c%.  with  a  small  amount  of  additive  noise  to 

account  for  the  residual  error.  The  usual  reason  for  transmitting  a 
"stretched  pulse"  instead  of  transmitting  p(t)  in  the  first  place  is  to 
maximize  signal  energy.  If  it  is  feasible  to  transmit  a  very  intense 
impulse  of  extremely  short  duration  then  matched  filtering  is 
unnecessary.  In  either  case,  the  sonar  output  is  thus  approximately  of 
the  form  of  Eq.  (19).  and  the  Vandermonde  factorization  method  can  be 
used  to  estimate  the  target  resonance  parameters. 

Although  the  form  of  Eq.  (19)  appears  different  from  that  for  which 
the  Vandermonde  factorization  was  applied,  it  is  actually  the 
same  as  Eq.  (1)  with  s^s-o^*  j  2Ttfk,  sk+£/2  =  sk* ,  dk  =  Akexp(j<pk),  and 

dk*K/2=d|(*  fork  =  1,2 . K/2.  A  target  with  K/2  resonances  thus  has  K 

spectral  components  or  "poles",  which  appear  in  conjugate  pairs.  In 
accordance  with  the  procedure  described  in  Section  Il.B  above,  we 
formed  the  data  matrix  A  whose  (m,n)th  element  is  defined  to  be 
a(m*n-2).  In  so  doing,  we  were  assuming  that  the  time  axis  has  been 
scaled  so  that  it  measures  t  in  units  of  the  data  sampling  interval.  In 
the  examples  to  be  described  below  we  assumed  a  sampling  interval  of 
1  ms;  i.e.,  t  expresses  time  in  units  of  milliseconds.  The  sampled  values 
of  a(t)  can  be  expressed  as 

K 

a(n)=  I  dkzkn  .  (21) 

k=1 


where 


zk  =  exp(sk)  =  exp(-o<k-j2TCfk)  .  (22) 

The  data  matrix  A  must  admit  the  Type  2  factorization  A  =  VC, 
where  the  zk's  appear  in  the  second  row  of  the  M*K  Vandermonde  V.  The 


target  resonance  parameters  cty  and  are  thus  determined  entirely  by 
|zk|  and  arg(Z|<),  and  the  accuracy  with  which  one  has  determined  them 
can  be  measured  in  terms  of  the  relative  error  in  |  |  and  the  angular 

error  in  arg(zj<). 

To  have  some  basis  of  comparison  with  existing  techniques  the  z^'s 
were  also  determined  by  a  least-squares  version  of  the  Prony  method; 
specifically,  the  equation  AATx  =  0  was  solved  for  an  Mxi  vector  x  (in 

practice,  by  finding  the  unit  vector  that  minimized  the  norm  of  AATX; 
see  Ref.  5)  and  the  z^'s  were  then  identified  as  a  subset  of  the  M-1 

roots  of  the  polynomial  (1,z,z2 . zM-1)x.  In  selecting  the  proper  roots 

it  was.  of  course,  helpful  to  know  the  correct  answers  (a  circumstance 
that  was  only  possible  because  the  data  were  artificially  generated  with 
known  resonance  parameters). 

Several  different  waveforms  parameter  sets  were  used  for  study 
purposes,  using  even  values  of  K  up  to  6.  To  summarize  the  results  we 
shall  discuss  the  performance  in  selected  cases  that  used  a  target 
transient  response  h(t)  comprising  two  strongly  damped  resonances 
(K=4)  at  f^  =  430  Hz  and  f2  =  280  Hz,  with  decay  times  of  (o^)-1  =  2  ms 
and  (o^)"1  =  3.333  ms,  and  with  initial  amplitudes  of  1.0  and  1.33, 
respectively. 

To  test  the  computational  methods  a  9x12  data  matrix  A  was  first 
constructed  from  h(t)  directly,  with  no  distortive  convolution.  As 
expected,  the  Vandermonde  factorization  method  identified  the  four  z^'s 
(in  two  conjugate  pairs)  with  great  precision.  The  least-squares  Prony 
method  produced  eight  z^'s,  of  which  half  were  at  the  correct  locations. 
Then  white  noise  was  added  to  the  sampled  data,  in  an  amount  adjusted 
to  perturb  the  computed  z^’s  slightly.  The  perturbations  in  the  z^'s 

were  virtually  the  same  for  both  the  Vandermonde  factorization  and 
least-squares  Prony  methods.  (The  observed  perturbations  in  the  higher 
resonance  frequency  were  proportionately  larger  than  those  of  the  lower 
resonance  frequency.) 


These  results  were  in  agreement  with  the  proffered  performance 
advantage  of  the  Vandermonde  factorization  method:  It  determines  the 
Z|<'s  with  a  degree  of  precision  comparable  to  that  obtained  with  a 
least-squares  Prony  analysis  of  higher  order,  but  without  the 
requirement  for  selecting  the  correct  z^'s  from  among  the  extraneous 
“noise  poles". 

Following  this  initial  test,  the  performance  was  tested  for 
convolved  versions  of  h(t),  using  a  15x36  form  of  the  data  matrix  A. 
The  results  for  three  special  cases  are  summarized  below.  In  each  case, 
convolution  was  performed  by  sampling  h(t)  and  p(t)  at  a  much  higher 
frequency  (10  kHz)  to  avoid  finite-sampling-rate  effects,  and  using 
409.6  ms  segments  of  data  to  avoid  windowing  effects.  Convolution 
was  actually  performed  by  using  an  8192-point  FFT. 

Case  1  It  was  assumed  that  the  sonar  transmitted  a  linear  FM 
“chirp'  of  61.5  ms  duration,  sweeping  from  150  Hz  to  800  Hz 
(timexbandwidth  =  40).  The  autocorrelation  function  of  this  chirp  (i.e., 
the  result  of  convolving  the  chirp  with  its  time-reversed  self)  became 
p(t);  its  approximate  duration  was  1.5  ms  (since  the  matched-filter- 
compressed  pulse  always  has  a  timexbandwidth  product  of  approximately 
one).  This  p(t)  was  convolved  with  h(t).  To  avoid  contaminating  the 
data  with  the  echo  components  still  being  generated  while  the  effective 
transmitted  pulse  p(t)  was  “active*,  the  initiation  of  data  sampling  was 
delayed  for  1.5  ms  beyond  the  onset  of  signal. 

The  Vandermonde  factorization  method  identified  the  higher 
resonance  frequency  well  (within  2*)  but  overestimated  its  decay  time 
constant  by  over  30*.  On  the  other  hand,  the  decay  time  constant  of  the 
lower  resonance  was  accurately  estimated,  but  its  frequency  was 
underestimated  by  about  15*. 


19 


The  seven  conjugate  pairs  of  z^’s  produced  by  the  least-squares 

Prony  method  were  scattered  around  the  left  half  of  the  unit  circle  in 
the  complex  z-plane;  however,  those  seen  to  be  closest  to  the  actual, 
known  values  were  closer  than  the  z^'s  measured  by  the  Vandermonde 

factorization  method. 

Case  2  For  the  next  exercise  it  was  assumed  that  the  sonar  system 
“forgot"  to  perform  the  matched  filtering  operation,  so  that  a(t)  was. 
itself,  a  61.5  ms  chirp  with  two  very  broad  resonance  peaks  (broad 
because  of  the  high  degree  of  damping)  at  the  two  resonance 
frequencies.  Thus  a(t)  did  not  even  remotely  approximate  Eq.  (19),  being 
dominated  by  the  FM  character  of  the  transmitted  pulse  rather  than  by 
the  system  resonances,  and  one  would  expect  that  the  z^'s  would  be 

improperly  identified.  Indeed,  the  least-squares  Prony  method  gave  a 
looping  pattern  of  z^'s  (as  viewed  in  the  complex  plane)  that  bore  no 
relation  to  the  actual  system  resonances.  Surprisingly,  however,  the 
Vandermonde  factorization  method  identified  both  resonance  frequencies 
(within  35K).  although  it  failed  to  identify  the  decay  rates  (indeed,  it 
perceived  the  resonances  to  be  undamped). 

The  ability  of  the  Vandermonde  factorization  method  to  find  the 
resonances  under  these  extraordinary  conditions  is  noteworthy,  but 
unexplained.  The  performance  continued  to  be  unexpectedly  good 
(although  not  always  as  good)  even  when  the  relationship  between 
sampling  rate  and  resonance  frequencies  was  varied,  and  the  number  of 
resonances  was  increased  to  three. 

Case  3  To  test  the  ability  of  the  Vandermonde  factorization  to 
accommodate  sharp-cutoff  bandpass  filtering  it  was  assumed  that  p(t) 
was  the  impulse  response  of  an  ideal  bandpass  filter  whose  cutoff 
frequencies  were  150  Hz  and  800  Hz,  the  same  as  the  limits  of  the  FM 
chirp  used  before.  Since  the  energy  spectrum  of  the  chirp  is  almost  flat 
over  this  frequency  range  and  small  elsewhere,  its  autocorrelation 
function  (which  is  the  inverse  Fourier  transform  of  the  energy  spectrum) 
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is  very  nearly  the  same  as  the  impulse  response  of  an  ideal  bandpass 
filter  (which  is  the  inverse  Fourier  transform  of  a  perfectly  flat,  band 
limited,  frequency  response;  i.e.,  it  is  a  modulated  "sine"  function). 
Therefore,  it  was  expected  that  the  results  would  be  very  similar  to 
those  of  Case  1 . 

However,  the  Vandermonde  factorization  failed  entirely  to  "see"  the 
higher  resonance  frequency,  and  estimated  both  of  its  resonance 
frequencies  to  be  near  the  lower  of  the  two  actual  resonance 
frequencies;  furthermore,  both  estimates  were  too  highly  damped  (one 
was  very  highly  damped).  On  the  other  hand,  the  least-square  Prony 
method  did  spectacularly  well  in  the  sense  that  a  subset  of  its  z^'s  fell 

almost  exactly  upon  the  correct  values. 

Perhaps  the  following  interpretation  explains  this  unexpected 
behavior:  Apparently  there  were  enough  extraneous  poles  (there  were 
10)  to  model  the  characteristics  of  the  bandpass  filter,  and  to  take  the 
"stress'  off  the  target  poles.  On  the  other  hand,  since  it  did  not  have 
these  extra  degrees  of  freedom,  the  Vandermonde  factorization  was 
forced  to  compromise  between  modeling  the  target  resonances  and 
modeling  the  bandpass  filter’s  characteristics,  with  disruptive  results. 

The  examples  of  transient  spectral  analysis  demonstrated  that 
Vandermonde  factorization  could  identify  spectral  components  of 
transients  of  very  short  duration.  However,  as  regards  analysis  of  sonar 
target  resonances,  the  scope  of  the  project  did  not  permit  an  extensive 
investigation  of  performance.  Only  a  preliminary  study  of  target 
resonance  analysis  could  be  made,  and  it  raised  as  many  questions  as  it 
answered.  The  surprising  ability  of  the  Vandermonde  factorization 
method  to  extract  target  resonances  buried  in  an  FM  chirp  of  long 
duration  was  even  more  curious  than  its  poor  performance  in  the 
presence  of  a  sharp-cutoff  bandpass  filter. 


As  the  broadband  direction  finding  application  was  studied  it 
appeared  more  and  more  likely  that  the  most  useful  applications  were 
those  of  small  order,  in  particular,  M=2.  For  that  case  only  one  sonar 
target  can  be  located  (i.e.,  K=1),  and  the  broadband  direction  finder  has 
little  to  offer  in  the  classical  case  of  passive  tracking  of  stationary 
targets  emitting  stationary,  broadband  random  process  of  very  low 
signal-to-noise  ratio.  For  that  environment  a  crosscorrelation  receiver 
pair  works  well,  and  has  the  ability  to  resolve  multiple  targets. 

However,  for  targets  emitting  brief,  albeit  intense,  transient  sounds 
the  broadband  direction  finder  is  useful,  since  it  can  be  configured  for 
M=2  if  it  can  be  assumed  that  the  target  emissions  are  so  brief  that 
only  one  occurs  at  a  given  time  (i.e.,  K=1  so  that  M=2  gives  M>K  as 
required). 

For  M=2  the  Vandermonde  Type  2  factorization  is  trivial.  (The 
singular  vectors  of  A  are  the  eigenvectors  of  the  2*2  matrix  AAH.  and 
therefore  a  simple  closed  form  solution  for  the  2*1  "matrix“  X  of 
Step  1.  All  that  is  required  to  get  P  is  to  scale  X  so  that  its  top 
element  is  1.  The  “polynomial"  of  Eq.  (13)  is  therefore  of  first  order.  It 
follows  that  a  closed  form  solution  for  z ^  (=u,)  can  be  written  in  terms 

of  the  three  distinct  elements  of  the  Hermitian,  2x2  matrix  AAH.  The 
algebra  is  straightforward,  if  tedious. 

However,  for  M=2  the  structure  is  simple  enough  that  an  alternate 
method  for  determining  the  direction  of  acoustic  transients  then 
becomes  attractive.  It  is  based  upon  the  observation  that  if  the  two 
outputs  b<|(t)  and  b2(t)  are  applied  to  the  horizontal  and  vertical  plates 
of  a  cathode  ray  tube  (CRT),  then  the  following  behavior  will  be 
exhibited:  In  the  absence  of  directional  signals  the  spot  will  hover  near 
the  center  of  the  screen,  producing  a  'fuzzy  spot"  due  to  random  noise. 
If  an  acoustic  transient  arrives  from  a  target  at  direction  cosine  u,  then 


the  displayed  spot  will  deflect  along  the  line  b2  =  ub1  in  accordance  with 
Eq.  (9).  The  length  of  the  trace  will  be  proportional  to  the  incident 
wave’s  amplitude,  but  the  incident  wave’s  spectrum  has  little  relevance; 
the  slope  u  indicates  target  direction. 

Suppose  that  the  observed  slope  is  “quantized’’;  i.e.,  the  CRT  is 
divided  into  a  large  number  of  uniformly  spaced  wedges,  each  covering  a 
small  span  of  values  of  the  slope  angle,  tan-1(u).  Suppose  that  for  each 
such  wedge  there  is  a  separate  power  detector  that  produces  no  output 
unless  the  spot  falls  into  the  appropriate  wedge,  whereupon  the  detected 
power  output  is  b12  +  b22.  If  all  of  these  zonal  power  detectors  are 
connected  through  a  bank  of  identical  smoothing  filters  to  a 
multichannel  recorder  it  will  provide  a  kind  of  histogram  record  of  the 
directional  energies  of  incident  wave  transients,  as  a  function  of  time. 

It  is  even  possible  to  scale  and/or  transform  the  display  coordinates 
to  “circularize"  the  background  noise’s  fuzzy  spot  so  that  the 
distribution  of  energies  due  to  background  noise  will  be  equal  in  all  of 
the  zones.  With  this  normalization,  the  presence  of  a  weak  target  can 
be  more  easily  perceived.  This  coordinate  transformation  has  a  mildly 
distortive  effect  upon  the  mapping  from  the  wedge  zones  to  actual 
target  angle,  but  is  not  a  serious  problem. 

In  a  rather  extensive  study,  this  concept  was  developed  in  detail  in 
a  full  technical  paper  that  was  submitted  for  publication  (see  Ref.  3). 
The  study  also  included  tests  of  performance  using  computer-generated 
data. 

The  results  were  very  encouraging,  and  indicated  that  this  proposed 
signal  processing  system  might  provide  an  attractive  alternative  for 
interception  of  directional  transients  of  arbitrary  spectra  and  brief 
duration,  with  enough  signal  energy  to  give  an  instantaneous  signal-to- 
noise  level  of  sufficient  size. 


IV.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  STUDY 


The  basic  principles  of  Vandermonde  factorization  were  developed 
without  encountering  any  fundamental  flaws  in  the  concept.  A  step-by- 
step  factorization  method  was  defined  and  tested  successfully.  The 
method  worked  when  applied  to  determining  target  resonances  in  an 
active  sonar  application  (using  computer-generated  data),  but  with  some 
unexpected  and  unexplained  behavior.  Our  limited  efforts  did  not  prove 
that  Vandermonde  factorization  is  the  method  of  choice  for  extracting 
target  resonances  from  echoes  elicited  by  sonar  transmissions  of  very 
brief  (either  before  or  after  pulse  compression)  duration. 

The  ability  of  the  Vandermonde  method  to  extract  target  resonances 
from  uncompressed  echoes  elicited  by  FM  chirp  transmissions  of  long 
duration  were  encouraging.  Moreover,  significant  progress  was  made  on 
the  problem  of  observing  brief  acoustic  transients  of  arbitrary  spectra, 
using  systems  of  such  low  order  that  the  actual  Vandermonde 
factorization  became  unnecessary,  so  that  a  directional-histogram 
approach  was  a  feasible  alternative.  The  results  of  this  new  approach 
to  acoustic  transient  interception  were  very  interesting. 

Recommended  for  further  research  are  the  following. 

(1)  Analysis  of  the  effects  of  noisy  data  upon  the  Vandermonde 
factorization. 

(2)  Development  of  a  well  documented  Vandermonde  factorization 
software  package  for  general  distribution. 

(3)  Further  studies  of  performance  of  Vandermonde  factorization  in 
target  resonance  analysis  for  both  compressed  pulses  and  uncompressed 
FM  chirps. 

(4)  Advanced  testing  of  the  directional-histogram  approach  to 
acoustic  transient  interception,  using  real  acoustic  data. 


REFERENCES 


T.  L.  Henderson,  “Rank  Reduction  for  Broadband  Waves  incident  on  a 
Linear  Receiving  Aperture,"  Applied  Research  Laboratories  Technical 
Paper  No.  85-21  (ARL-TP-85-21 ).  Applied  Research  Laboratories, 
The  University  of  Texas  at  Austin,  1985. 

T.  L.  Henderson,  R.  S.  Bailey,  and  S.  G.  Lacker,  "Use  of  Rectangular 
Vandermonde  Matrix  Decomposition  To  Resolve  Damped  Sinusoidal 
Oscillation  Components."  J.  Acoust.  Soc.  Am.  80(1 )  (1986). 

T.  L.  Henderson.  "A  Broadband  Bearing  Deviation  Indicator  for 
Acoustic  Transients,"  Applied  Research  Laboratories  Technical  Paper 
No.  86-20  (ARL-TP-86-20),  Applied  Research  Laboratories,  The 
University  of  Texas  at  Austin.  1986. 

S.  Haykin  (ed.),  Methods  of  Nonlinear  Spectral  Analusis.  2nd  edition, 
(Springer- Verlag,  New  York,  1983). 

T.  L.  Henderson,  “Geometric  Methods  for  Determining  System  Poles 
from  Transient  Response,"  IEEE  Trans.  Acoustics,  Speech,  and  Signal 
Processing,  ASSP-29.  No.  5  (1981). 

G.  W.  Stewart,  Introduction  to  Matrix  Computation  (Academic  Press, 
New  York,  1973). 


DISTRIBUTION  LIST  FOR 
ARL-TR-87-29 

UNDER  CONTRACT  N0001 4-85-K-0501 

Codu  No. 

Office  of  Naval  Research 
Department  of  the  Navy 
Arlington,  VA  2221 7 

1  -  2  Attn:  R.  Fitzgerald  (Code  1 125AO) 

Director 

Naval  Research  Laboratory 
Washington,  DC  20375 
3  Attn:  Code  2627 

4-15  Commanding  Officer  and  Director 

Defense  Technical  Information  Center 
Bldg.  5,  Cameron  Station 
5010  Duke  Street 
Alexandria,  VA  22314 

16  Advanced  Sonar  Division,  ARL:UT 

1 7  C.  Robert  Culbertson,  ARL:UT 

18  Terry  L.  Henderson.  ARL:UT 

19  Library,  ARL:UT 


20  -  30 


Reserve,  ARL:UT 


